function [ll, grad] = T1FAVParametersLL_m( par, x, echoes, sigma, scales)

% computes the likelihood of the rice distributed Flip Angle Variation data
% echoes = angles [n x 1] vector
% par = the 2 element parametervector [M,E]
%
% The datamodel is given by
% SI(i) = M(1-E)*sin(echoes(i)) / (1-E*cos(echoes(i)))
% dSI/dM = (1-E)*sin(alp)/(1-E*cos(alp))
%
%          (M*sin)(E*cos-1)-(E*sin-1)*M*cos
% dSI/dE = ---------------------------------
%                  (1-E*cos)^2
%
% dM/dPar = {dSI/dM, dSI/dE}
%
% Created by Henk Smit, EMC, 01-2011 based on the work by Dirk Poot, University of Antwerp, 13-8-2007

if size(par,1) ~= 2% + (nargin<=3)  HENK
    error('incorrect parameter vector');
end;
numtr = size(par,2);
if nargin<=3
    sigma = par(3,:); %HENK
end;
numgr = size(echoes,1); %numAngles
if ((size(x,1)~=numgr || size(x,2)~=numtr) && ~isempty(x)) %|| size(echoes,2)~=1  || numel(sigma)~=1
    error('incorrect size in input.');
end;

SI = zeros(numgr,numtr);
dSIdM = SI;
dSIdE = SI;

%scales=[1000;1000;0.001]; %stefan: to bring the parameters in the same range
%scales=[1;1;1];
%par=par./scales;

for k=1:numtr
    M=par(1,k);
    E=par(2,k);
    angles=echoes(:,1);
    dSIdM(:,k) = ((1-E).*sin(angles))./(1-E.*cos(angles));
    dSIdE(:,k) = ((sin(angles).*cos(angles).*M.*(1-E)) - (M.*sin(angles).*(1-E.*cos(angles))))./ ( (1 - (E.*cos(angles))).^2);
    SI = M*dSIdM;
end;

if isempty(x)
    ll = SI;
    return;
end;

if nargout>1 %henk scaling
    [lrpdf, ricegrad] = logricepdf(x, SI, sigma,logical([0 1 nargin<=3])); 
    dAdpar=[dSIdM dSIdE];
    grad = reshape(-sum( reshape( ricegrad(:,ones(1,2)).*dAdpar, numgr, numtr * 2) ), numtr, 2);
    grad = grad';
   
    %grad = grad';
else
    [lrpdf] = logricepdf(x, SI, sigma );
end;
ll = - sum( lrpdf(:) );
